{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 一、梯度下降（Gradient Descent）小结\n",
    "\n",
    "1. 求解模型参数时，即无约束优化问题时，常采用梯度下降和最小二乘。\n",
    "2. 梯度向量，从几何意义上讲，就是函数变化增加最快的地方（方向）。\n",
    "3. 代数方法描述、矩阵方法描述 、算法调优-归一化、步长、初始值\n",
    "4. 批量梯度下降、随机梯度下降、小批量梯度下降（Mini-batch）\n",
    "5. 其他无约束优化算法，最小二乘（求解析解，需要求逆矩阵）、牛顿法和拟牛顿法   （二阶的海森矩阵的逆矩阵或伪逆矩阵求解） \n",
    "\n",
    "## 二、算法描述\n",
    "\n",
    "### 2.1 前提假设\n",
    "\n",
    "1. 步长（Learning rate）：梯度前进的长度\n",
    "2. 特征（feature）与样本  ：$(x^{(0)},y^{(0)}),(x^{(1)},y^{(1)})$   \n",
    "3. 假设函数（hypothesis function）：$h_{\\theta}(x) = \\theta_0+\\theta_1x$\n",
    "\n",
    "4. 损失函数（loss function） ：模型拟合好坏， 对于m个样本$（x_i,y_i）(i=1,2,...m)$，采用线性回归均方差，损失函数为：$J(\\theta_0, \\theta_1) = \\sum\\limits_{i=1}^{m}(h_\\theta(x_i) - y_i)^2$\n",
    "\n",
    " 　其中$x_i$表示第$i$个样本特征，$y_i$表示第$i$个样本对应的输出，$h_θ(x_i)$为假设函数。 \n",
    "\n",
    "### 2.2 算法过程\n",
    "\n",
    "初始化$\\theta$参数,步长$\\alpha$\n",
    "\n",
    "1. 确定损失函数的梯度$\\frac{\\partial}{\\partial\\mathbf\\theta}J(\\mathbf\\theta)$\n",
    "2. 步长乘以梯度，得到当前位置下降的距离 $\\alpha\\frac{\\partial}{\\partial\\theta}J(\\theta)$\n",
    "\n",
    "3. 更新$\\theta$,更新完后判断截止条件，回到步骤一。\n",
    "   $$\n",
    "   \\mathbf\\theta= \\mathbf\\theta - \\alpha\\frac{\\partial}{\\partial\\theta}J(\\mathbf\\theta)\n",
    "   $$\n",
    "\n",
    "### 2.3 矩阵法描述\n",
    "\n",
    "假设函数：$h_\\mathbf{\\theta}(\\mathbf{X}) = \\mathbf{X\\theta}$\n",
    "\n",
    "损失函数：$J(\\mathbf\\theta) = \\frac{1}{2}(\\mathbf{X\\theta} - \\mathbf{Y})^T(\\mathbf{X\\theta} - \\mathbf{Y})$\n",
    "\n",
    "$h_\\mathbf{\\theta}(\\mathbf{X})、Y$：mx1维      $X$：mx(n+1)维     $\\theta$：(n+1)x1维       $J(\\mathbf\\theta)$：标量  1x1维\n",
    "\n",
    "\n",
    "\n",
    "1. 损失函数梯度：$\\frac{\\partial}{\\partial\\mathbf\\theta}J(\\mathbf\\theta) = \\mathbf{X}^T(\\mathbf{X\\theta} - \\mathbf{Y})$  利用链式法则和复合函数求导（也可以利用维度相容原则）\n",
    "\n",
    "2. 步长乘以梯度，得到当前位置下降的距离 $\\alpha\\frac{\\partial}{\\partial\\theta}J(\\theta)$\n",
    "\n",
    "3. 更新$\\theta$,更新完后判断截止条件，回到步骤一。\n",
    "   $$\n",
    "   \\mathbf\\theta= \\mathbf\\theta - \\alpha\\frac{\\partial}{\\partial\\theta}J(\\mathbf\\theta)\n",
    "   $$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 三、代码实现验证"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "元素乘法：np.multiply(a,b) 、  *   \n",
    "矩阵乘法：np.dot(a,b) 或 np.matmul(a,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "optimal: [[0.51583286]\n",
      " [0.96992163]]\n",
      "error function: [[1.01496241]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1b52c66e688>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3AAAAHSCAYAAACtoSkbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzde5jXY+L/8efdVIwc0pZDOS87WIdisNZ3WXbXOC1hlxxjkdiwDoNS67edmUJCckoIOSSsXeO4DotlElsOI4fQDB1kEkam6f798Z5sMqlpDu/PZ+b5uK65Zj735/Rq5pJec9/v+w4xRiRJkiRJma9V2gEkSZIkSavGAidJkiRJWcICJ0mSJElZwgInSZIkSVnCAidJkiRJWcICJ0mSJElZonXaAWrTsWPHuMUWW6QdQ5IkSZJSMWXKlHkxxk7Lj2dkgdtiiy0oKSlJO4YkSZIkpSKE8GFt4y6hlCRJkqQsYYGTJEmSpCxhgZMkSZKkLGGBkyRJkqQsYYGTJEmSpCxhgZMkSZKkLGGBkyRJkqQsYYGTJEmSpCxhgZMkSZKkLGGBkyRJkqQsYYGTJEmSpCxhgZMkSZKkLGGBkyRJkqQsYYGTJEmSpCxhgZMkSZKkLNE67QCSJEmS1JQmTy2jqLiU8opKOrfPpbAgj+7duqQda5VY4CRJkiS1GJOnltF30jQqq6oBKKuopO+kaQBZUeJcQilJkiSpxSgqLv2uvC1VWVVNUXFpSonqxgInSZIkqcUor6is03imscBJkiRJajE6t8+t03imscBJkiRJajEKC/LIbZPzvbHcNjkUFuSllKhu3MREkiRJUouxdKMSd6GUJEmSpCzQvVuXrClsy3MJpSRJkiRlCQucJEmSJGUJC5wkSZIkZQkLnCRJkiRlCQucJEmSJGUJC5wkSZIkZQkLnCRJkiRlCQucJEmSJGUJC5wkSZIkZQkLnCRJkiRlCQucJEmSJGWJlRa4EMKmIYSnQwhvhRDeCCGcUzPeIYTweAhhRs3n9Vfw/J41j5kRQujZ0H8ASZIkSWopVmUGbjFwfoxxO+AXwJ9DCNsDFwNPxhi3AZ6suf09IYQOwKXAHsDuwKUrKnqSJEmSpB+30gIXY/wkxvhqzdcLgbeALsBhwPiah40Hutfy9ALg8Rjj/Bjj58DjwAENEVySJEmSWpo6XQMXQtgC6Ab8B9gwxvgJJCUP2KCWp3QBPl7m9qyaMUmSJElSHa1ygQshrA3cD/wlxvjFqj6tlrG4gtfvFUIoCSGUzJ07d1VjSZIkSdLq+fLLtBPU2SoVuBBCG5LyNiHGOKlmeHYIYeOa+zcG5tTy1FnApsvc3gQor+09Yow3xBjzY4z5nTp1WtX8kiRJklQ3H34IZ54JXbrAJ5+knaZOVmUXygDcDLwVY7ximbseApbuKtkTeLCWpxcD+4cQ1q/ZvGT/mjFJkiRJalrvvw+nnQZbbw033QTHHAOhtkWDmav1KjxmL+AEYFoI4bWasX7AcOCeEMIpwEfAHwFCCPlA7xjjqTHG+SGEQcArNc8bGGOc36B/AkmSJEn6Me+8A0OHwh13QOvWcPrpcNFFsOmmK39uhgkx1npJWqry8/NjSUlJ2jEkSZIkZbM334QhQ+Duu2GNNZLiVlgInTunnWylQghTYoz5y4+vygycJEmSJGWPadNg8GC4915Yay04//zkY8MN005WbxY4SZIkSc3D1KkwaBA88ACssw707QvnngsdO6adrMFY4CRJkiRlt5dfTorb3/8O660Hl14KZ58NHTqknazBWeAkSZIkZacXXoCBA6G4OClrgwdDnz5JiWumLHCSJEmSssszzyTF7amnoFMnGD48OddtnXXSTtboLHCSJEmSMl+M8OSTyVLJZ59NNiQZOTLZWbJdu7TTNRkLnCRJkqTMFSM8+mhS3F58Ebp0gauvhlNPhdzctNM1uVZpB5AkSZKkH4gRHn4Ydt8dDjoIyspgzBh47z0466wWWd7AAidJkiQpkyxZApMmwS67wKGHwmefwY03wowZ0Lt3ciB3C2aBkyRJkpS+6mqYOBF23hmOPBK++gpuvRVKS5Plkm3bpp0wI1jgJEmSJKVn8WK44w7YYQfo0SMpchMmwFtvQc+e0KZN2gkzigVOkiRJUtOrqoJx42C77eCEE5Kids89MH06HHss5OSknTAjuQulJEmSpAY1eWoZRcWllFdU0rl9LoUFeXTv1iW589tvYfx4GDoUZs6Ebt2Sa94OOwxaOb+0MhY4SZIkSQ1m8tQy+k6aRmVVNQBlFZX0nTSNVou+4dBXi5NDtz/+ONldcvRoOPhgCCHl1NnDAidJkiSpwRQVl35X3gDWrPqGHiXF7Hn1JPjiM/jlL5NdJfff3+K2GixwkiRJkhpMeUUlAGt9W8lxU/9Jr1cm0emrCl7adAc6PXAP7Luvxa0eLHCSJEmSGszWuZHfPn0fp778AD+p/ILnNu/Knw/rQdmOu/Hv/fZLO17Ws8BJkiRJqr+KChg9mkeuGEnbLxbw9Fa7MvqXPXi1y3bktslhWEFe2gmbBQucJEmSpNU3fz5cdRWMGgVffEHbQw/lX3/sRf9ZuZRXVNJl+V0oVS8WOEmSJEl1N3cuXHEFXHMNfPklHHkk9O8PXbvya+DfaedrpixwkiRJklbdp5/CiBEwZgxUVsJRRyXFbYcd0k7WIljgJEmSJK1cWRkUFcHYsclh3MceC/36wXbbpZ2sRbHASZIkSVqxjz6Cyy6Dm26C6mo48UTo2xe22SbtZC2SBU6SJEnSD33wAQwbBrfemtw++WS4+GLYcstUY7V0FjhJkiRJ/zNjBgwdCrffDjk50KsXXHghbLZZ2smEBU6SJEkSwFtvwZAhcNdd0LYt9OkDhYXQxe3/M4kFTpIkSWrJpk+HwYPhnnsgNxfOOw/OPx822ijtZKqFBU6SJElqiV57DQYNgkmTYO21k+vbzj0XOnVKO5l+hAVOkiRJakleeSUpbg8/DOutB3/9K5xzDnTokHYyrQILnCRJktQSvPgiDBwIjz4K66+flLg+faB9+7STqQ4scJIkSVJz9uyzSVl74gno2DE5GuDMM2HdddNOptVggZMkSZKamxjh6aeTGbdnnoENN4QRI6B3b2jXLu10qgcLnCRJktRcxAiPPZYUtxdegM6dYdQoOO20ZIdJZT0LnCRJkrScyVPLKCoupbyiks7tcyksyKN7tww+Dy1GeOSRpLi98gpsuilcdx2cfDKsuWba6dSALHCSJEnSMiZPLaPvpGlUVlUDUFZRSd9J0wAyr8QtWQIPPphc4zZ1Kmy5JdxwA/TsmRzGrWanVdoBJEmSpExSVFz6XXlbqrKqmqLi0pQS1aK6Ojl4u2tXOOIIWLgQxo2D0tJkuaTlrdmywEmSJEnLKK+orNN4k6quhjvvhB13hKOPhqoquOMOeOstOOkkaNMm7YRqZBY4SZIkaRmd29e+2ceKxptEVRWMHw/bbQfHHQc5OTBxIkyfntxu7ZVRLYUFTpIkSVpGYUEeuW1yvjeW2yaHwoK8pg/z7bdw002Ql5fMsLVrB/ffD6+/DkcdlRQ5tShWdUmSJGkZSzcqSXUXykWL4JZbYPhw+Ogj2G235DiAQw6BEJouhzKOBU6SJElaTvduXdLZcbKyEm68ES6/HMrKYM89YexYKCiwuAmwwEmSJEnp++oruP56KCqC2bNh772Ta97228/ipu9ZaYELIdwCHALMiTHuUDM2EVi6CLg9UBFj7FrLc2cCC4FqYHGMMb+BckuSJEnZb+HC5MDtESNg3jz4zW+SzUn22SftZMpQqzIDdytwDXDb0oEY49FLvw4hjAQW/Mjz940xzlvdgJIkSVKzs2ABjB4NV14J8+fDAQfAgAHwy1+mnUwZbqUFLsb4bAhhi9ruCyEE4Chgv4aNJUmSJDVD8+cnm5GMGpWUuN//Hvr3h913TzuZskR9r4H7FTA7xjhjBfdH4LEQQgTGxhhvqOf7SZIkSdln3jy44gq45ppk2eThhyfFbZdd0k6mLFPfAncMcNeP3L9XjLE8hLAB8HgI4e0Y47O1PTCE0AvoBbDZZpvVM5YkSZKUAWbPTq5vGzMGvv4a/vjHpLjtuGPayZSlVvsg7xBCa+AIYOKKHhNjLK/5PAd4AFjh3HCM8YYYY36MMb9Tp06rG0uSJElKX3k5nHsubLllMvPWvTu88UayQYnlTfWw2gUO+C3wdoxxVm13hhDahRDWWfo1sD8wvR7vJ0mSJGW2jz+GPn1gq62STUqOPhrefhvuuAO22y7tdGoGVlrgQgh3AS8CeSGEWSGEU2ru6sFyyydDCJ1DCP+oubkh8HwI4XXgZeCRGOOjDRddkiRJyhAzZ8Lpp8NPf5ocvH3CCfDOOzBuHGyzTdrp1Iysyi6Ux6xg/KRaxsqBg2q+fh/YuZ75JEmSpMz17rswdCjcfju0agWnngoXXQSbb552MjVT9d3ERJIkSWp5SkthyBCYMAHatoUzz4QLL4QuXdJOpmbOAidJkiStqjfegMGDk81IcnOTjUouuAA22ijtZGohLHCSJElqcpOnllFUXEp5RSWd2+dSWJBH924ZPHv12mtJcbv/flh77WS27bzzYIMNUomTdd8/NRgLnCRJkprU5Kll9J00jcqqagDKKirpO2kaQOaVkJISGDQIHnoI1l03OcPtL3+Bn/wktUhZ9f1Tg6vPMQKSJElSnRUVl35XPpaqrKqmqLg0pUS1eOklOPhg2G03ePZZ+Nvf4MMPkzKXYnmDLPn+qdE4AydJkqQmVV5RWafxJvXcc0lJe/zxpKgNGZKc67buumkn+05Gf//U6JyBkyRJUpPq3D63TuONLkZ4+mnYd1/Ye294/XUoKkrOduvXL6PKG2Tg909NygInSZKkJlVYkEdum5zvjeW2yaGwIK9pg8QIjz0Gv/oV7LdfcjTAlVfCBx8kO0uuvXbT5llFGfP9UypcQilJkqQmtXSjjdR2UYwR/vEPGDgQXn4ZNtkErrkGTjkF1lyzaTLUQ+rfP6UqxBjTzvAD+fn5saSkJO0YkiRJak6WLEl2kxw0CF59FTbfPFki2bMnrLFG2umk7wkhTIkx5i8/7hJKSZIkNW9LlsC990K3bnD44bBgAdxyC8yYAb16Wd6UVSxwkiRJap6qq+Guu2DHHeGoo2DRIrjtNnj7bTj5ZGjTJu2EUp1Z4CRJktS8LF6cFLXtt4djj03G7roL3ngDTjgBWrsNhLKXBU6SJEnNw7ffws03Q15ecl1bbi7cdx9MmwY9ekBOzspfQ8pw/vpBkiRJ2W3RIrj1Vhg2DD78EHbdFR58EH7/ewgh7XRSg3IGTpIkSdnpm2+S7f+33hp694aNNoJHHoFXXoFDD7W8qVlyBk6SJEnZ5euvYexYuPxy+PRT+L//S3aV/O1vLW1q9ixwkiRJyg5ffgnXXQcjR8KcObDvvsnmJPvsY3FTi2GBkyRJUmb74otkqeQVV8Bnn8H++8OAAcnMm9TCWOAkSZKUmT7/HK6+Gq66Cioq4OCDk+K2xx5pJ5NSY4GTJElSZvnsM7jyShg9Opl9O+ywpLjtumvayaTUWeAkSZKUGebMSa5vu+46+OorOPJI6N8fdt457WRSxrDASZIkNUOTp5ZRVFxKeUUlndvnUliQR/duXdKOVbtPPoERI2DMmORogB494JJL4Oc/TzuZlHEscJIkSc3M5Kll9J00jcqqagDKKirpO2kaQGaVuFmzkqMAbrgBFi+G446Dfv0gLy/tZFLG8iBvSZKkZqaouPS78rZUZVU1RcWlKSVazocfwhlnwE9/msy6HXccvP02jB9veZNWwhk4SZKkZqa8orJO403mvfdg2LCkqIUAf/oTXHwxbLFFurmkLGKBkyRJamY6t8+lrJay1rl9bgppgHfegSFDYMIEaN0aeveGCy+ETTdNJ4+UxVxCKUmS1MwUFuSR2ybne2O5bXIoLGji5YlvvgnHHgvbbQf33gtnnw0ffJAcD2B5k1aLM3CSJEnNzNKNSlLbhfK//4XBg+G++2CtteD885OPDTdsmveXmjELnCRJUjPUvVuXpt9x8tVXYdAgmDwZ1lkH+vaFc8+Fjh2bNofUjFngJEmSVD//+U9S3B55BNq3h0svhXPOgfXXTzuZ1OxY4CRJkrR6/v1vGDgQHnsMOnRIlk326QPrrZd2MqnZssBJkiRp1cUIzzyTFLenn4ZOneCyy5Jz3dZZJ+10UrNngZMkSdLKxQhPPJEslXzuOdhoIxg5Ek4/Hdq1Szud1GJY4CRJkrRiMcKjjyYzbi+9BF26wNVXw6mnQm5K58pJLZjnwEmSJOmHYoSHHoLdd4eDDoLychgzBt57D846y/ImpcQCJ0mSpP9ZsgTuvx+6dYPDDoP58+Gmm2DGDOjdG9ZYI+2EUotmgZMkSRJUV8Pdd8NOO8Ef/gCVlTB+PJSWwimnQNu2aSeUhAVOkiSpZVu8GG6/HX7+czjmmGQG7s474c034cQTobVbJkiZxAInSZLUElVVwS23wLbbJkVtjTXgnntg+vSkyOXkpJ1QUi38lYokSVJLsmhRsjRy2DCYORN22QUeeAAOPRRa+bt9KdP5X6kkSVJL8M03cO21sPXWydltG24If/87lJRA9+6WNylLrPS/1BDCLSGEOSGE6cuM/b8QQlkI4bWaj4NW8NwDQgilIYR3QwgXN2RwSZIkrYKvv4arroKttoI+fWDzzaG4GF58EQ4+GEJIO6GkOliVX7XcChxQy/iVMcauNR//WP7OEEIOcC1wILA9cEwIYfv6hJUkSdIq+vJLKCqCLbeEc8+FvDx46il47jnYf3+Lm5SlVlrgYozPAvNX47V3B96NMb4fY/wWuBs4bDVeR5IkSavqiy+S69u22AIuvDA5FuDZZ+Hpp2HffS1uUparz2LnPiGE/9YssVy/lvu7AB8vc3tWzVitQgi9QgglIYSSuXPn1iOWJElSC1RRAQMHJsWtXz/YY49kmeTjj8OvfpV2OkkNZHUL3Bjgp0BX4BNgZC2Pqe3XO3FFLxhjvCHGmB9jzO/UqdNqxpIkSWphPvsMBgxIrm279FLYe2945RV45BH4xS/STiepga3WMQIxxtlLvw4h3Aj8vZaHzQI2Xeb2JkD56ryfJEmSljN3Lowcmews+eWXcOSR0L8/dO2adjJJjWi1ClwIYeMY4yc1Nw8HptfysFeAbUIIWwJlQA/g2NVKKUmSpMSnn8KIETBmDFRWwtFHwyWXwA47pJ1MUhNYaYELIdwF/BroGEKYBVwK/DqE0JVkSeRM4PSax3YGbooxHhRjXBxC6AMUAznALTHGNxrlTyFJktTclZXB5ZfDDTfAt9/Ccccl17ptu23aySQ1oRDjCi9LS01+fn4sKSlJO4YkSVL6PvoIhg+Hm2+GJUvgxBOhb9/kQG5JzVYIYUqMMX/58dVaQilJkqRG9v77yXEA48cnt08+GS6+ODnXTVKLZYGTJEnKJO+8A0OHwh13QOvW0KsXXHQRbLrpyp8rqdmzwEmSJGWCt96CIUPgrrugbVs46ywoLITOndNOJimDWOAkSZLSNG0aDB4M994Lublw3nlw/vmw0UZpJ5OUgSxwkiRJaZg6FQYNggcegHXWSa5vO/dc6NQp7WSSMpgFTpIkqSm98kpS3B5+GNZbD/76VzjnHOjQIe1kkrKABU6SJKkpvPBCUtwefRTWXz/5uk8faN8+7WSSsogFTpIkqTE980xS1p58Ejp2TM50O/PMZNmkJNWRBU6SJKmhxQhPPQUDB8Kzz8KGG8KIEdC7N7Rrl3Y6SVnMAidJktRQYoTi4mTG7YUXkiMARo2C005LdpiUpHpqlXYASZKkrBdjsinJHnvAgQfCxx/DddfBe+/B2Wdb3iQ1GAucJEnS6lqyJDkGYNdd4dBDYd48uPFGePddOOMMWHPNtBNKamYscJIkSXVVXQ333AM77wxHHAELF8K4cVBaCqeeCm3bpp1QUjNlgZMkSVpVixfDhAmwww5w9NHJ7dtvh7fegpNOgjZt0k4oqZmzwEmSJK1MVRXceitsvz0cfzy0bg0TJ8L06f+7LUlNwL9tJEmSVuTbb+G222DoUPjgA+jaFSZNgsMOg1b+HlxS0/NvHkmSpOUtWgRjxsA22yRHAHTsCA89BK++CocfbnmTlBpn4CRJkpaqrEx2kbzsMigvhz33hLFjoaAAQkg7nSRZ4CRJkvjqK7j+eigqgtmzYZ99kqWT++1ncZOUUSxwkiSp5Vq4MDlwe8SI5Ay33/wm2Zxkn33STiZJtbLASZKklmfBAhg9Gq68EubPhwMOgAED4Je/TDuZJP0oC5wkSWo55s+HUaOSjwUL4Pe/h/79Yffd004mSavEAidJkpq/efPgiivgmmuSZZNHHJEUt27d0k4mSXVigZMkSc3X7NnJ9W1jxsDXX8NRR8Ell8COO6adTJJWiwVOkiQ1P+XlyY6SY8cmZ7odc0xS3LbbLu1kklQvFjhJktR8fPxxcobbTTfB4sVwwgnQr19yILckNQMWOEmSlP1mzoRhw2DcuOT2SSfBxRfDVlulmUqSGpwFTpIkZa93302K2223QatWcNppcNFFsNlmaSeTpEZhgZMkSdmntBSGDIEJE6BtW/jzn6GwELp0STuZJDUqC5wkScoeb7wBgwfDxImQmwvnngsXXAAbbZR2MklqEhY4SZKU+V5/HQYNgvvvh7XXTpZJnncedOqUdjJJalIWOEmSlLmmTEmK24MPwrrrwoABcM458JOfpJ1MklJhgZMkSZnnpZf4tPASNnr+KSrWXJt7f3sSnQcUcvDe26ed7DuTp5ZRVFxKeUUlndvnUliQR/duXoMnqXFZ4CRJUuZ47rlkxu3xx1kjd10u26cnt3c7mC/XWIvcxz+iap31MqIkTZ5aRt9J06isqgagrKKSvpOmAWREPknNV6u0A0iSpBYuRnj6adh3X9h7b3j9da458HT26n0zY37xR75cYy0AKquqKSouTTlsoqi49LvytlQm5ZPUfFngJElSOmKExx5LStt++yVHA1x1FXzwASN3+j1ft839wVPKKypTCPpDK8qRKfkkNV8WOEmS1LRihEcegT33hIIC+PBDuPZaeP/9ZIOStdaic/sfljdgheNNLdPzSWq+LHCSJKlpLFkCkydDfj4ccgjMng033ADvvgtnnglrrvndQwsL8shtk/O9p+e2yaGwIK+pU9cq0/NJar7cxESSJDWuJUuS89sGD4b//he23hpuuQWOPx7atKn1KUs3AsnUXR4zPZ+k5ivEGNPO8AP5+fmxpKQk7RiSJKk+qqvhnnuS4vbmm7DttnDJJdCjB7T2d8iS9GNCCFNijPnLj7uEUpIkNazFi+G222D77eHYYyEEuPtumD49mXWzvEnSaltpgQsh3BJCmBNCmL7MWFEI4e0Qwn9DCA+EENqv4LkzQwjTQgivhRCcUpMkqTn79lu4+WbIy4OePSE3F+67L1k2efTRkJOz8teQJP2oVZmBuxU4YLmxx4EdYow7Ae8AfX/k+fvGGLvWNv0nSZKagUWL4Prr4Wc/g1NPhQ4d4MEHYepUOPJIaOWCH0lqKCv9GzXG+Cwwf7mxx2KMi2tuvgRs0gjZJElSJqushNGj4ac/hTPOgI03hn/8A15+GQ49NFk6KUlqUA3xK7E/Af9cwX0ReCyEMCWE0KsB3kuSJKXtq6/giitgq63g7LOTz48/Di+8AAceaHGTpEZUr6uIQwiXAIuBCSt4yF4xxvIQwgbA4yGEt2tm9Gp7rV5AL4DNNtusPrEkSVJjWLgQxoyBESNg7lzYb79kc5J99kk7mSS1GKs9AxdC6AkcAhwXV3AWQYyxvObzHOABYPcVvV6M8YYYY36MMb9Tp06rG0uSJDW0BQtgyBDYYgu46CLYZRd4/nl48knLmyQ1sdWagQshHABcBOwTY/x6BY9pB7SKMS6s+Xp/YOBqJ5UkSU3r889h1Kjko6ICDjkEBgyA3Vf4+1hJUiNbaYELIdwF/BroGEKYBVxKsuvkGiTLIgFeijH2DiF0Bm6KMR4EbAg8UHN/a+DOGOOjjfKnkCRJDWfePLjyymSDkoULoXt36N8fdt017WSS1OKttMDFGI+pZfjmFTy2HDio5uv3gZ3rlU6SJDWd2bNh5Ei47jr4+mv4wx+S4rbTTmknkyTVqNcmJpIkqRn45BMoKkrOclu0CHr0gEsuge23TzuZJGk5FjhJklqqWbPgssvgxhth8WI4/njo1y85kFuSlJEscJIkNZLJU8soKi6lvKKSzu1zKSzIo3u3LmnHgpkzYfhwGDcOliyBnj2hb9/kQG5JUkazwEmS1AgmTy2j76RpVFZVA1BWUUnfSdMA0itx770HQ4fCbbdBq1bwpz/BxRfD5punk0eSVGerfQ6cJElasaLi0u/K21KVVdUUFZc2fZjS0mSWLS8PJkyAM85IytyYMZY3ScoyzsBJktQIyisq6zTeKN54IzmAe+JEWGMNOOccuOAC2HjjpssgSWpQzsBJktQIOrfPrdN4g3r9dfjjH2HHHeGhh5LSNnNmckSA5U2SspoFTpKkRlBYkEdum5zvjeW2yaGwIK/x3nTKlOTQ7a5dobg42VFy5sxkp8kNNmi895UkNRmXUEqS1AiWblTSJLtQ/uc/MGgQPPIItG8P/+//wdlnw/rrN/x7SZJSZYGTJKmRdO/WpXF3nHz++aS4PfYYdOiQXO/25z/Deus13ntKklJlgZMkKZvECM88AwMHwtNPQ6dOyRLJM86AddZJO50kqZFZ4CRJygYxwhNPJDNuzz0HG20EV1wBvXpBu3Zpp5MkNRELnCRJmSxG+Oc/k+L20kvQpQuMHg2nnAK5TbCjpSQpo7gLpSRJmSjG5AiA3XaDgw+G8nK4/vrkAO4+fSxvktRCWeAkScokS5bA/fdDt25w2GHw+edw000wYwacfnpyILckqcWywEmSlAmqq+Huu2GnneAPf4DKShg/HkpLk+WSbdumnVCSlAEscJIkpWnxYrjjDvj5z+GYY5IZuDvvhDffhBNPhNZeri5J+h8LnCRJaaiqgnHjYNtt4YQTkqWR99wD06cnRS4nJ+2EkqQM5K/1JElqSosWJUs5anEAACAASURBVEsjhw2DmTNhl13ggQfg0EOhlb9XlST9OP9PIUlSU/jmG7j2Wth662Qzkg02gL//HUpKoHt3y5skaZU4AydJUmP6+mu48Ua47DL45BPYay+4+Wb43e8ghLTTSZKyjAVOkqTG8OWXybltRUUwZw78+tcwYULy2eImSVpNFjhJ0gpNnlpGUXEp5RWVdG6fS2FBHt27dUk7Vmb74otkqeTIkfDZZ8lM24AB8KtfpZ3sB/z5SlL2scBJkmo1eWoZfSdNo7KqGoCyikr6TpoG4D/ya1NRAVdfDVddlRy+fdBBSXH7xS/STlYrf76SlJ28YlqSVKui4tLv/nG/VGVVNUXFpSklylCffZYUtc03h0svTWbaXnkFHnkkY8sb+POVpGzlDJwkqVblFZV1Gm9x5s5Nlklee21yvduRR0L//tC1a9rJVok/X0nKTs7ASZJq1bl9bp3GW4xPP4ULLoAttoDLL4dDDoFp0+C++7KmvIE/X0nKVhY4SVKtCgvyyG2T872x3DY5FBbkpZQoZWVlcM45sOWWcOWVyYzbm2/CXXfBDjukna7O/PlKUnZyCaUkqVZLN7Jo8bsUfvQRDB+enN22ZAmccAL065ccyJ3F/PlKUnYKMca0M/xAfn5+LCkpSTuGJKkle/99GDYMxo9Pbp98Mlx8cTIDJ0lSIwshTIkx5i8/7gycJEnLmjEDhg6F22+H1q2hVy+46CLYdNO0k0mSZIGTJAmAt96CIUOSa9rWWAPOOgsKC6Fz57STSZL0HQucJKllmzYNBg+Ge++F3Fw4//zkY8MN004mSdIPWOAkSS3T1KkwaBA88ACss05yfdu550KnTmknkyRphSxwkqSW5ZVXkuL28MOw3nrw178mxwN06JB2MkmSVsoCJ0lqGV54ISlujz6alLVBg5Lr3NZbL+1kkiStMgucJKl5e+aZpKw9+SR07Jic6XbmmcmySUmSsowFTpLU/MQITz0FAwfCs88mG5KMHAmnnw7t2qWdTpKk1WaBkyQ1HzFCcXEy4/bCC8kRAKNGwWmnJTtMSpKU5VqlHUCSpHqLMdmUZI894MADYdYsuO46eO89OPtsy5skqdmwwEmSsteSJTBpEuy6Kxx6KMybBzfeCDNmwBlnwJprpp1QkqQGZYGTJGWf6mq45x7YeWc48kj48ku49VYoLYVTT4W2bdNOKElSo1ilAhdCuCWEMCeEMH2ZsQ4hhMdDCDNqPq+/guf2rHnMjBBCz4YKLklqgRYvhgkTYIcd4OijkyI3YQK8+Sb07Alt2qSdUJKkRrWqM3C3AgcsN3Yx8GSMcRvgyZrb3xNC6ABcCuwB7A5cuqKiJ0nSClVVJTNs228Pxx+fFLWJE2HaNDj2WGjtnlySpJZhlQpcjPFZYP5yw4cB42u+Hg90r+WpBcDjMcb5McbPgcf5YRGUJKl2336bXNOWlwcnnwxrr51c8/baa3DUUZCTk3ZCSZKaVH2ugdswxvgJQM3nDWp5TBfg42Vuz6oZkyRpxb75JtlFcuutoVev5ADuhx+GKVPg8MOhlZdwS5Japsb+P2CoZSzW+sAQeoUQSkIIJXPnzm3kWJKkjFRZCVdfDT/9Kfz5z7DppvDoo/Cf/8Ahh0Co7X8rkiS1HPUpcLNDCBsD1HyeU8tjZgGbLnN7E6C8theLMd4QY8yPMeZ36tSpHrEkSVnnq69g5EjYcks45xzYZht48kl4/nkoKLC4SZJUoz4F7iFg6a6SPYEHa3lMMbB/CGH9ms1L9q8ZkyQJFi6E4cNhiy3gggtgxx3hmWfgX/+C/fazuEmStJxVPUbgLuBFIC+EMCuEcAowHPhdCGEG8Lua24QQ8kMINwHEGOcDg4BXaj4G1oxJklqyigoYNCgpbn37Qn4+/Pvf8PjjsPfeaaeTJCljhRhrvSQtVfn5+bGkpCTtGJKkhjZ/PowalXwsWAC//z0MGAC77ZZ2MkmSMkoIYUqMMX/5cQ/OkSQ1vnnz4Ior4JprkmWTRxwB/ftDt25pJ5MkKatY4CRJjWf2bBgxIjkSoLIyObvtkkuSa90kSVKdWeAkSQ2vvBwuvxzGjk0O4z7mmKS4bbdd2skkScpqFjhJUsP56CO47DK4+WZYvBhOOAH69UuOBZAkSfVmgZMk1d8HHyTHAYwbl9w+6SS4+GLYaqtUY0mS1NxY4CRJq+/dd2HoULjtNsjJgdNOg4sugs02SzuZJEnNkgVOklR3b78NQ4bAnXdC27bQpw8UFkKXLmknkySpWbPASZJW3fTpMHgw3HMP5ObCeefB+efDRhulnUySpBbBAidJWrnXXoNBg2DSJFh77WSZ5HnnQadOaSeTJKlFscBJklaspCQpbg89BOuuCwMGwDnnwE9+knYySZJaJAucJOmHXnwxKW7//Cesvz787W9w9tnQvn3aySRJatEscJKk/3nuORg4EJ54Ajp2hGHD4Mwzk9k3SZKUOgucJLV0McLTTyczbv/6F2ywARQVQe/eyfVukiQpY1jgJKmlihEeeywpbv/+N2y8MVx1VXKW21prpZ1OkiTVolXaASRJTSxGeOQR+MUv4IAD4KOP4Npr4f33kw1KLG+SJGUsC5wktRRLlsDkyZCfD4ccAnPmwA03wLvvJte5rblm2gklSdJKWOAkqblbsgTuvRe6doXDD4cvvoBx4+Cdd5Llkm3bpp1QkiStIgucJDVX1dVw552w445w1FFQVQW33w5vvQUnnQRt2qSdUJIk1ZEFTpKam8WL4bbbYPvt4bjjIAS4+26YPh2OPx5au3+VJEnZygInSc3Ft9/CzTdDXh707Am5uXD//fDf/8LRR0NOTtoJJUlSPVngJCnbLVoE118P22wDp54KHTrAQw/B1KlwxBHQyr/qJUlqLlxHI0nZqrISbroJLrsMysqSYwHGjoWCgmTZpCRJanYscJKUbb76KilqRUXw6afwq1/BrbfCb35jcZMkqZmzwElStli4EMaMgREjYO5c2G+/ZHOSffZJO5kkSWoiFjhJynQLFsA118AVV8D8+ckSyQEDYK+90k4mSZKamAVOkjLV55/DqFHJR0UFHHJIUtx23z3tZJIkKSUWOEnKNPPmwZVXwujRybLJww+H/v1hl13STiZJklJmgZOkTDFnDowcCddeC19/DX/8I1xyCey0U9rJJElShrDASVLaPvkk2VHy+uuTM9169EiK2/bbp51MkiRlGAucJKVl1qzkDLcbb4TFi+H446FfP/jZz9JOJkmSMpQFTpKa2syZMHw4jBsHS5bASSfBxRfDT3+adjJJkpThLHCS1FTeew+GDYPx46FVKzjlFLjoIth887STSZKkLGGBk6TGVloKQ4fChAnQpg2ceSYUFsImmzB5ahlFdz1FeUUlndvnUliQR/duXdJOnDUmTy2jqLjU758kqcWwwElSY3njDRgyBCZOhDXWgHPOgQsugI03BpLy0XfSNCqrqgEoq6ik76RpAJaQVeD3T5LUErVKO4AkNTuvv54cAbDjjvDQQ8ls28yZyREBNeUNoKi49LvysVRlVTVFxaVNHDg7+f2TJLVEzsBJUkOZMgUGDYIHH4R1102OAvjLX+AnP6n14eUVlXUa1/f5/ZMktUTOwElSff3nP3DIIZCfD888A3/7G3z4YVLmVlDeADq3z63TuL7P758kqSWywEnS6nr+eSgogF/8Al56Kdmo5MMP4a9/hfbtV/r0woI8ctvkfG8st00OhQV5jZW4WfH7J0lqiVxCKUl1EWMyyzZwIDz9NGywAVx+OZxxBqy9dp1eaulGG+6iuHr8/kmSWqIQY0w7ww/k5+fHkpKStGNI0v/ECE88kRS3559PNiO58ELo1QvWWivtdJIkqZkJIUyJMeYvP+4SSkn6MTHCP/4Be+4J+++f7CZ5zTXw/vvJBiWWN0mS1IQscJJUmxiT3SR32w0OPhg+/RSuvx7efRf+/GdYc820E0qSpBZotQtcCCEvhPDaMh9fhBD+stxjfh1CWLDMY/5a/8iS1IiWLIH77oNu3aB7d6iogJtvhhkz4PTTkwO5JUmSUrLam5jEGEuBrgAhhBygDHigloc+F2M8ZHXfR5KaRHU13HMPDBkCb7wBeXlw221wzDHQ2v2eJElSZmioJZS/Ad6LMX7YQK8nSU1j8WK4/Xb4+c/h2GOTpZN33ZWUuBNOsLxJkqSM0lAFrgdw1wru2zOE8HoI4Z8hhJ+v6AVCCL1CCCUhhJK5c+c2UCxJWoGqKrjlFth2WzjxxGRp5L33wrRp0KMH5OSs/DUkSZKaWL0LXAihLXAocG8td78KbB5j3BkYDUxe0evEGG+IMebHGPM7depU31iSVLtFi2DsWNhmGzjllOTA7cmTYepU+MMfoJV7O0mSpMzVEP9SORB4NcY4e/k7YoxfxBi/rPn6H0CbEELHBnhPSaqbb75Jtv/femvo3Rs22ggeeQReeQUOO8ziJkmSskJDXNxxDCtYPhlC2AiYHWOMIYTdSQrjZw3wnpK0ar7+OplxKyqCTz6B//u/ZOnkb38LIaSdTpIkqU7qVeBCCGsBvwNOX2asN0CM8XrgD8AZIYTFQCXQI8YY6/OekrRKvvwSxoyBESNgzhz49a/hzjthn30sbpIkKWvVq8DFGL8GfrLc2PXLfH0NcE193kOS6uSLL5KlkldcAZ99BvvvDwMGJDNvkiRJWc79sSU1D59/DldfDVddlRy+ffDBSXHbY4+0k0mSJDUYC5yk7PbZZ3DllTB6dDL7dthhSXHbdVcAJk8to6i4lPKKSjq3z6WwII/u3bqkHFqSJGn1WOAkZac5c2DkSLjuOvjqKzjySOjfH3be+buHTJ5aRt9J06isqgagrKKSvpOmAVjiJElSVnLfbEnZ5ZNP4PzzYYstkp0lf//75PDte+/9XnkDKCou/a68LVVZVU1RcWkTBpYkSWo4zsBJyg6zZsHll8ONN0JVFRx7LPTrB9tuu8KnlFdU1mlckiQp01ngJGW2Dz+E4cOTs9uWLIETT4S+fZMDuVeic/tcymopa53b5zZGUkmSpEbnEkpJmen99+G005KidvPNcPLJMGNG8vUqlDeAwoI8ctvkfG8st00OhQV5jZFYkiSp0TkDJymzvPMODB0Kd9wBrVtD795w4YWw6aZ1fqmlG5W4C6UkSWouLHCSMsObb8KQIXD33bDGGnD22XDBBdC5c71etnu3LhY2SZLUbFjgJKXrv/+FwYPhvvtgrbWSHSbPPx823DDtZJIkSRnHAicpHVOnwqBB8MADsM46ycYk554LHTumnUySJCljWeAkNa2XX06K29//DuutB5demiyX7NAh7WSSJEkZzwInqWm88AIMHAjFxUlZGzwY+vRJSpwkSZJWiQVOUuN65pmkuD31FHTqlJzpduaZybJJSZIk1YkFTlLDixGefDJZKvnss8mGJCNHwumnQ7t2aaeTJEnKWhY4SQ0nxmSJ5MCB8OKLyREAV18Np54Kublpp5MkScp6rdIOIKkZiBEefhh23x0OPBDKyuC66+C99+CssyxvkiRJDcQCJ2n1LVkCkybBLrvAoYfCZ5/BjTfCjBlwxhmw5pppJ5QkSWpWLHCS6q66GiZOhJ13hiOPhK++gnHjoLQ0WS7Ztm3aCSVJkpolC5ykVbd4MUyYADvsAD16JEVuwgR480046SRo0ybthJIkSc2aBU7SylVVwa23wnbbwfHHJ0Vt4kSYNg2OPRZaux+SJElSU/BfXZJW7NtvYfx4GDoUZs6Ebt2Sa94OOwxa+fsfSZKkpua/wCT90DffJLtIbr019OoFG2yQ7DI5ZQocfrjlTZIkKSXOwEn6n8pKuOEGuPxyKC+HX/4y2VVy//0hhLTTSZIktXgWOEnJLpLXXw9FRTB7NuyzD9x+O+y7r8VNkiQpg1jgpJZs4UK49loYORLmzYPf/hbuuQf23jvtZJIkSaqFBU5qiSoqYPRouOoqmD8fDjwQBgyAPfdMO5kkSZJ+hAVOaknmz4dRo5KPBQvg0EOhf3/Ybbe0k0mSJGkVWOCklmDePLjiCrjmmmTZ5BFHJMWtW7e0k0mSJKkOLHBSczZ7NowYkRwJUFkJRx0Fl1wCO+6YdjJJkiStBguc1ByVlydHAYwdmxzGfeyx0K8fbLdd2skkSZJUDxY4qTn56CO47DK4+WZYvBhOPBH69oVttkk7mSRJkhqABU5qDj74AIYPh3HjktsnnwwXXwxbbpluLkmSJDUoC5yUzd59F4YOhdtug5wcOO00uOgi2GyztJNJkiSpEVjgpGz09tswZAjceSe0bQt9+kBhIXTpknYySZIkNSILnJRNpk9PitvEiZCbC+edB+efDxttlHYySZIkNQELnJQNXnsNBg+G+++HtddOrm8791zo1CntZJIkSWpCFjgpk5WUwKBB8NBDsN56MGAA/OUv0KFD2skkSZKUAguclIlefDEpbv/8J6y/PgwcCGedBe3bp51MkiRJKbLASZnkueeSsvbEE9CxIwwbBmeeCeuum3YySZIkZQALnJS2GOHpp5MZt3/9CzbcEEaMgN69oV27tNNJkiQpg9S7wIUQZgILgWpgcYwxf7n7AzAKOAj4Gjgpxvhqfd9XynoxwmOPJcXt3/+Gzp3hqquSs9zWWivtdJIkScpADTUDt2+Mcd4K7jsQ2KbmYw9gTM1ntRCTp5ZRVFxKeUUlndvnUliQR/dumXNeWZPnixEeeSQpbi+/DJtuCtdeC3/6E6y5ZuO9ryRJkrJeUyyhPAy4LcYYgZdCCO1DCBvHGD9pgvdWyiZPLaPvpGlUVlUDUFZRSd9J0wAyosQ1ab4lS5LdJAcOhKlTYcst4YYboGfP5DBuSZIkaSVaNcBrROCxEMKUEEKvWu7vAny8zO1ZNWNqAYqKS78rR0tVVlVTVFyaUqLva5J8S5bAvfdC165w+OGwcCGMGwelpclyScubJEmSVlFDzMDtFWMsDyFsADweQng7xvjsMveHWp4Tlx+oKX+9ADbbbLMGiKVMUF5RWafxptao+aqrYeLE5ADut96CbbeFO+6Ao4+G1u4fJEmSpLqr9wxcjLG85vMc4AFg9+UeMgvYdJnbmwDltbzODTHG/BhjfqdOneobSxmic/vcOo03tUbJt3gxjB8P228Pxx0HOTlw990wfXpy2/ImSZKk1VSvAhdCaBdCWGfp18D+wPTlHvYQcGJI/AJY4PVvLUdhQR65bXK+N5bbJofCgryUEn1fg+b79lu46SbIy4OTTkp2krz/fnj99WTWLSdnpS8hSZIk/Zj6TgVsCDyQnBRAa+DOGOOjIYTeADHG64F/kBwh8C7JMQIn1/M9lUWWbgSSqbtQNki+RYuSa9qGDYOPPoL8/OQ4gEMOgVDbCmJJkiRp9YRkc8jMkp+fH0tKStKOIf24yspkxu2yy6CsDPbcE/76VygosLhJkiSpXkIIU5Y/Yxua5hgBqXn56isYOxaKiuDTT2HvvZNr3vbbz+ImSZKkRmWBk1bVwoVw3XUwciTMnQu/+U2yOck++6SdTJIkSS2EBU5amQULYPRouPJKmD8/WSI5YADstVfaySRJktTCWOCkFZk/H0aNSj4WLEg2JRkwAHZf/qQMSZIkqWlY4KTlzZuXzLaNHp0smzz8cOjfH3bZJe1kkiRJauEscNJSs2cn17dddx18/TX88Y9wySWw005pJ5MkSZIAC5wE5eXJjpJjxyZnuvXokRS37bdPO5kkSZL0PRY4tVwff5yc4XbTTbB4MRx/PPTrBz/7WdrJJEmSpFpZ4NTyzJwJw4fDLbdAjHDSSdC3L2y1VdrJJEmSpB9lgVPL8e67MGwY3HYbtGoFp54KF10Em2+edjJJkiRplVjg1PyVlsKQITBhArRtC2eeCYWFsMkmaSeTJEmS6sQCp+brjTdg8GCYOBHWXBP+8he44ALYeOO0k0mSJEmrxQLXDEyeWkZRcSnlFZV0bp9LYUEe3bt1STtWel5/PSlu990H7drBhf+/vTsPsqss8zj+fWgC1bIYkc0giyVTwcoIBCMuCBUYIBqRZAacCliIwgihBmQzSBJAJUQgYRPQYCZQbBFlCZESMyHWoKNYUGZBgyQxIsHJIkGSGFNpNAnP/HFuqKbtm62777m37/dTleruc87tfrrf+55zfnnfc84VcNllsO++ZVfWK/n+kyRJqh0DXIObPm8Zo6fNp23DJgCWrWlj9LT5AM13Ej1nDowbBz/8Iey5Z/Hw7UsugXe/u+zKei3ff5IkSbW1U9kFqGsmzlz01snzZm0bNjFx5qKSKirBs8/Cpz8NgwbBz34G3/gGvPJKEeYMbz3K958kSVJtOQLX4Javaduu5b3KL34B114Ls2YVQW38eLjwwmL0TTXR1O8/SZKkEjgC1+D69W3druUNLxOefhqOPx6OPba43m3ChOLZbmPGGN5qrOnef5IkSSUzwDW4UUP609qn5W3LWvu0MGpI/5Iq6iGZxUjbccfBCSfAwoVw663w8svFIwF2373sCptS07z/JEmS6oRTKBvc5htF9Nq7AGbCjBnFVMnnniue3XbHHXDuudDqKE/Zev37T5Ikqc5EZpZdwz8YNGhQzp49u+wyVKZMeOKJ4kYkc+bAwQfD6NHwhS/ArruWXZ0kSZLUoyJiTmYO6rjcKZSqL2++WTy/beBAGD4cVq+Gu++GxYvh/PMNb5IkSWpqBjjVh02b4KGH4PDD4bOfhbY2uP9+WLQIzjkH+vQpu0JJkiSpdAY4lWvjRnjgARgwAM48s5g6+dBD8OKLcNZZsLOXaUqSJEmbGeBUjg0b4J574LDD4POfL6ZGPvIIzJ8PI0ZAS8vWv4ckSZLUZBzeUG397W9w771www3Fs9uOOgqmT4fPfAZ28v8TJEmSpC3xjFm18cYbcOedcOihMHIk7LcfPPkkzJ4Nw4YZ3iRJkqRt4Aicetb69TB5MkyYACtWwDHHFFMnTzwRIsquTpIkSWooBjj1jHXrYNIkuOkmWLkSBg+GBx+E4483uEmSJEk7yACn7rV2bTFV8pZb4PXX4aST4Oqr4dhjy65MkiRJangGOHWP1avh9tvhtttgzRoYOrQIbh/9aNmVSZIkSb2GAU5d8/rrcOutcMcdxejbsGFFcPvQh8quTJIkSep1DHDaMStXws03w3e+U1zvdvrpcNVVcMQRZVcmSZIk9VoGuG0wfd4yJs5cxPI1bfTr28qoIf0ZPvCAsssqx4oVxY1JJk0qHg0wYgSMHQsDBpRdmSRJktTrGeC2Yvq8ZYyeNp+2DZsAWLamjdHT5gM0V4hburR4FMDkybBxI5x5JowZA4cdVnZlkiRJUtPw6clbMXHmorfC22ZtGzYxceaikiqqsVdegQsugPe/vxh1+9znYOFCuP9+w5skSZJUY47AbcXyNW3btbzXeOkluP56uO++4rlt55wDV14JhxxSdmWSJElS0zLAbUW/vq0s6ySs9evbWkI1NfC738H48TB1Kuy8M4wcCVdcAQceWHZlkiRJUtNzCuVWjBrSn9Y+LW9b1tqnhVFD+pdUUQ958cXiurYPfAAeeQS+/GX4wx+KxwMY3iRJkqS64AjcVmy+UUmvvQvlb34D110Hjz4K73gHXH558W+//cquTJIkSVIHBrhtMHzgAb0nsG02dy6MGwfTp8Mee8Do0XDppbD33mVXJkmSJKkKA1yzee65Irg9+SS8853wta8V0yX32qvsyiRJkiRthQGuWTzzDFx7LTz1VBHWrrsOLrywCHGSJEmSGsIO38QkIg6MiKcjYkFE/DYiLu5km8ER8ZeIeL7y75qulavtkgk//SmccAJ84hMwbx7ceCMsWQJjxxreJEmSpAbTlRG4jcDlmTk3IvYA5kTErMx8scN2P8/MU7rwc7S9MuEnPymmSv7857D//nDzzXD++bDbbmVXJ0mSJGkH7fAIXGauyMy5lc//CiwAetmdPhpMJsyYAR//OJx8cvEYgNtvLz5edpnhTZIkSWpw3fIcuIg4BBgIPNfJ6o9FxK8jYkZEDNjC9zgvImZHxOzXXnutO8pqHpnwxBNw9NEwdCgsXw6TJsFLL8FFF0FrL33ouCRJktRkuhzgImJ34DHgksxc22H1XODgzDwCuAOYXu37ZObkzByUmYP22WefrpbVHN58Ex57DAYOhGHDYNUqmDIFFi+GkSNh113LrlCSJElSN+pSgIuIPhThbWpmTuu4PjPXZua6yuc/BvpEhA8a66pNm+D734fDD4fTT4f16+Hee2HhQjj3XNhll7IrlCRJktQDunIXygDuBhZk5i1Vttm/sh0RcXTl572+oz+z6W3cCA8+CAMGwBlnFCNwU6fCggVw9tnQp0/ZFUqSJEnqQV25C+UxwFnA/Ih4vrJsDHAQQGbeBZwOXBARG4E2YERmZhd+ZnPasKEIbuPHF9e1ffCD8PDDcNppsFO3XMYoSZIkqQHscIDLzF8AsZVt7gTu3NGf0fT+/vdiauT11xfPbhs4EB5/HE491eAmSZIkNSFTQD164w349rfh0EOLZ7ftuy/86EcwZw4MH254kyRJkppUV6ZQqrutXw+TJ8OECbBiRfE8tylT4KSTILY42ClJkiSpCRjg6sG6dXDXXTBxIqxcCYMHF9e8HX+8wU1Mn7eMiTMXsXxNG/36tjJqSH+GDzyg7LIkSZJUAgNcmdauLaZK3nIL/PnPcOKJcPXVcNxxZVemOjF93jJGT5tP24ZNACxb08boafMBDHGSJElNyIupyrBmDYwbB4ccAmPGwIc/DL/8JcyaZXjT20ycueit8LZZ24ZNTJy5qKSKJEmSVCZH4Gpp1Sq47Tb41reK0bdTTy1G3AYNKrsy1anla9q2a7kkSZJ6NwNcLbz2WjFN8s47i+vdTjsNrroKjjyy7MpU5/r1bWVZJ2GtX9/WEqqRJElS2ZxC2ZP+9Cf4yleKqZI33ginnALz58OjjxretE1GDelPa5+WrpFQaQAAB4BJREFUty1r7dPCqCH9S6pIkiRJZXIEricsW1bcUfK73y0exn3mmTB2LBx2WNmVqcFsvlGJd6GUJEkSGOC61x//WIy0TZkCb74JZ51V3KTk0EPLrkwNbPjAAwxskiRJAgxw3ePll+H66+Hee4uvv/hFuPJKeN/7Si1LkiRJUu9igOuKxYvhm9+EBx6AlhY47zy44go46KCyK5MkSZLUCxngdsTChTB+PHzve7DLLnDRRTBqFPTrV3ZlkiRJknoxA9z2eOEFuO46ePhhaG2Fyy6Dyy+H/fcvuzJJkiRJTcAAty1Wr4YvfQkeewx23724vu3SS2GffcquTJIkSVITMcBtiz33hOXL4Zpr4OKLYa+9yq5IkiRJUhMywG2LlhZ45hmIKLsSSZIkSU1sp7ILaBiGN0mSJEklM8BJkiRJUoMwwEmSJElSgzDASZIkSVKDMMBJkiRJUoMwwEmSJElSgzDASZIkSVKDMMBJkiRJUoMwwEmSJElSgzDASZIkSVKDMMBJkiRJUoMwwEmSJElSgzDASZIkSVKDMMBJkiRJUoMwwEmSJElSgzDASZIkSVKDMMBJkiRJUoMwwEmSJElSg4jMLLuGfxARrwGvlF2H2Bv4c9lF6C22R32xPeqL7VFfbI/6YnvUF9ujvtRzexycmft0XFiXAU71ISJmZ+agsutQwfaoL7ZHfbE96ovtUV9sj/pie9SXRmwPp1BKkiRJUoMwwEmSJElSgzDAaUsml12A3sb2qC+2R32xPeqL7VFfbI/6YnvUl4ZrD6+BkyRJkqQG4QicJEmSJDUIA1yTi4gDI+LpiFgQEb+NiIs72WZwRPwlIp6v/LumjFqbRUQsiYj5lb/17E7WR0TcHhG/j4jfRMRRZdTZDCKif7v3/fMRsTYiLumwjf2jB0XEPRGxMiJeaLdsr4iYFRGLKx/fVeW1Z1e2WRwRZ9eu6t6rSntMjIiFlf3R4xHRt8prt7hv0/ar0h5fj4hl7fZJQ6u89pMRsahyLLmydlX3XlXa4wft2mJJRDxf5bX2j25W7Ry3NxxDnELZ5CLiPcB7MnNuROwBzAGGZ+aL7bYZDHwlM08pqcymEhFLgEGZ2ekzSSoH44uAocBHgG9l5kdqV2FziogWYBnwkcx8pd3ywdg/ekxEHAesA+7PzH+uLJsArMrMGyonnu/KzK92eN1ewGxgEJAU+7YPZebqmv4CvUyV9jgZ+J/M3BgRNwJ0bI/KdkvYwr5N269Ke3wdWJeZN23hdS3A74CTgKXAr4Az2h/7tf06a48O628G/pKZ13aybgn2j25V7RwX+AINfgxxBK7JZeaKzJxb+fyvwALggHKr0lYMozg4ZGY+C/St7KTUs/4FeKl9eFPPy8z/BVZ1WDwMuK/y+X0UB+SOhgCzMnNV5YA7C/hkjxXaJDprj8x8KjM3Vr58FnhvzQtrUlX6x7Y4Gvh9Zv4hM/8OfJ+iX6kLttQeERHAvwMP1bSoJraFc9yGP4YY4PSWiDgEGAg818nqj0XEryNiRkQMqGlhzSeBpyJiTkSc18n6A4D/a/f1UgzdtTCC6gde+0dt7ZeZK6A4QAP7drKN/aQc5wAzqqzb2r5N3efCypTWe6pMD7N/1N6xwKuZubjKevtHD+pwjtvwxxADnACIiN2Bx4BLMnNth9VzgYMz8wjgDmB6retrMsdk5lHAp4D/rEzJaC86eY1zoXtQROwCnAo80slq+0d9sp/UWESMBTYCU6tssrV9m7rHJOD9wJHACuDmTraxf9TeGWx59M3+0UO2co5b9WWdLKubPmKAExHRh+KNPTUzp3Vcn5lrM3Nd5fMfA30iYu8al9k0MnN55eNK4HGKqS7tLQUObPf1e4HltamuaX0KmJuZr3ZcYf8oxaubpw1XPq7sZBv7SQ1VLvA/BfhcVrm4fhv2beoGmflqZm7KzDeB/6Lzv7P9o4YiYmfg34AfVNvG/tEzqpzjNvwxxADX5Cpzsu8GFmTmLVW22b+yHRFxNMX75vXaVdk8ImK3yoW2RMRuwMnACx02ewL4fBQ+SnFB9Ioal9psqv7Pqf2jFE8Am+8Idjbww062mQmcHBHvqkwhO7myTN0sIj4JfBU4NTPXV9lmW/Zt6gYdron+Vzr/O/8K+KeIeF9lhsEIin6lnnEisDAzl3a20v7RM7Zwjtvwx5Cdyy5ApTsGOAuY3+7WtmOAgwAy8y7gdOCCiNgItAEjqv0Pq7psP+DxSh7YGfheZv53RIyEt9rjxxR3oPw9sB74Ykm1NoWIeAfFndrOb7esfXvYP3pQRDwEDAb2joilwNeAG4CHI+Jc4I/AZyvbDgJGZuZ/ZOaqiBhHcaIKcG1m7sjNHtROlfYYDewKzKrsu57NzJER0Q+YkplDqbJvK+FX6FWqtMfgiDiSYrrXEir7rvbtUblj6IUUJ6QtwD2Z+dsSfoVepbP2yMy76eQaavtHTVQ7x234Y4iPEZAkSZKkBuEUSkmSJElqEAY4SZIkSWoQBjhJkiRJahAGOEmSJElqEAY4SZIkSWoQBjhJkiRJahAGOEmSJElqEAY4SZIkSWoQ/w80AnkYhl4VVAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "m = 20 #设置数据集大小\n",
    "X0 = np.ones((m,1))   #m行1列\n",
    "X1 = np.arange(1,m+1).reshape(m,1) #m行1列\n",
    "X = np.hstack((X0,X1)) #沿着水平方向将矩阵堆叠起来\n",
    "\n",
    "y = np.array([\n",
    "    3, 4, 5, 5, 2, 4, 7, 8, 11, 8, 12,\n",
    "    11, 13, 13, 16, 17, 18, 17, 19, 21\n",
    "]).reshape(m, 1)\n",
    "\n",
    "alpha = 0.01# 梯度\n",
    "\n",
    "\n",
    "#误差函数 均方误差 \n",
    "def error_function(theta,X,y):\n",
    "    diff = np.dot(X,theta)-y  \n",
    "    return (1./2/m)*np.dot(np.transpose(diff),diff) \n",
    "\n",
    "\n",
    "def gradient_function(theta,X,y):\n",
    "    diff = np.dot(X,theta)-y     #20*2与2*1相乘 = 20*1\n",
    "    return (1/m)*np.dot(np.transpose(X),diff)  #2*20与20*1 = 2*1\n",
    "\n",
    "def gradient_descent(X,y,alpha):\n",
    "    theta = np.array([1,1]).reshape(2,1)  #定义θ参数   起始参数自定义的为（1,1）\n",
    "    gradient = gradient_function(theta,X,y)\n",
    "    error = error_function(theta,X,y)\n",
    "    while not np.all(np.absolute(gradient) <= 1e-5):  #只要gradient的绝对值不小于1e-5   （一般判定条件是误差小于多少）\n",
    "        theta = theta - alpha*gradient\n",
    "        gradient = gradient_function(theta,X,y)\n",
    "    return theta       \n",
    "\n",
    "optimal = gradient_descent(X,y,alpha)\n",
    "print('optimal:',optimal)\n",
    "print('error function:',error_function(optimal,X,y)) \n",
    "\n",
    "fig,ax = plt.subplots(figsize=(15,8)) \n",
    "ax.scatter(X1,y)\n",
    "\n",
    "model_y = optimal[0]*X0 + optimal[1]*X1\n",
    "ax.plot(X1,model_y,color = 'red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "参考：[深入浅出--梯度下降法及其实现](https://www.jianshu.com/p/c7e642877b0e)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
